Counterflow instability and turbulence in a spin-1 spinor Bose— Einstein condensate 
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We theoretically study counterflow instability and turbulence in a spin-1 spinor Bose-Einstein 
condensate by the Gross-Pitaevskii equation and the Bogoliubov-de Gennes equation. Our study 
considers (i) the dynamics induced by the counterflow of two components with different magnetic 
quantum numbers, which leads to turbulence with spin degrees of freedom, and (ii) the properties 
of the turbulence. For (i), the behavior of the condensate induced by the counterflow strongly 
depends on whether the spin-dependent interaction is ferromagnetic or antiferromagnetic, leading 
to different behaviors for the dispersion relation and the spin density vector, etc. For (ii), we 
numerically calculate the spectrum of the spin-dependent interaction energy, which also depends 
on the spin-dependent interaction. The spectrum of the spin-dependent interaction energy in the 
ferromagnetic case exhibits a —7/3 power law, whereas that in the antiferromagnetic case does not. 
The —7/3 power law can be explained by scaling analysis. 

PACS numbers: 67.85.De,03.75.Lm,67.25.dk,47.37.+q 



I. INTRODUCTION 

Turbulence is one of the most important topics in mod- 
ern physics. Turbulence in classical fluids has been stud- 
ied for a long time p), but turbulence is also observed 
in diverse fields such as low- temperature physics, plasma 
physics, and astronomy, etc. Kolmogorov studied tur- 
bulence in classical fluids in 1941 [2] and significantly 
advanced the study of turbulence. He found that the ki- 
netic energy spectrum obeyed a —5/3 power law, which 
became known as the Kolmogorov —5/3 power law. 

Several methods have been used for generating turbu- 
lence in classical fluids. Reynolds performed a famous 
experiment in 1883 in which he created turbulence by 
using flow through a circular pipe. Another method uses 
hydrodynamic instabilities, which have been investigated 
for a long time. There are many kinds of instabilities such 
as Kelvin-Helmholtz (KH) and Rayleigh-Taylor (RT) in- 
stabilities [3| which can lead to turbulence. 

Hydrodynamic instability has recently been studied 
in atomic Bose-Einstein condensates (BECs). Atomic 
BECs exhibit different dynamics from classical fluids be- 
cause they are quantum fluids. Quantized vortices occur 
in BECs; these vortices are nucleated through the hydro- 
dynamic instability. KH and RT instabilities have been 
studied in two-component atomic BECs. For KH insta- 
bility, quantized vortices nucleate at the boundary layer 
between the two components [4], whereas RT instabil- 
ity generates the shape of mushroom in the condensate 
winded by the quantized vortex ring [5]. Moreover, coun- 
terflow instability in two-component atomic BECs results 
in vortices that nucleate by the collapse of solitons [6|, l7j . 
Thus, hydrodynamic instability has been extensively in- 
vestigated in two-component atomic BECs. 

Unlike two-component atomic BECs, spinor BECs 
have spin degrees of freedom and exhibit phenomena 
characteristic of spin. Collisions of spin-1 spinor BECs 



have been investigated numerically and show different 
results from one- and two-component BECs [8]. The 
hydrodynamic equation for spinor BECs has recently 
been studied [9l-[llj. However, hydrodynamic instability 
in spinor BECs has been investigated less than in two- 
component BECs. We expect that hydrodynamic insta- 
bility in spinor BECs will exhibit unique behavior due to 
their spin degrees of freedom and form turbulent states 
in which the spin density vector has different directions. 

Quantum turbulence has been investigated for a long 
time in superfluid 4 He and 3 He I12L being currently 
studied in atomic BECs too |13l-[l6j. Numerical stud- 
ies predict that the Kolmogorov —5/3 power law, which 
was first observed in classical fluids, also holds in atomic 
BECs [15| . Turbulence in two-component BECs has been 
investigated jQ, • We expect that the turbulent state in 
a spin-1 spinor BEC will exhibit properties characteris- 
tic of the spin degrees of freedom, which one- and two- 
component BECs do not show. This is one of the major 
themes of this paper. 

In this paper, we focus on counterflow instability in 
spin-1 spinor BECs in a homogeneous two-dimensional 
system. There are three reasons for studying the insta- 
bility in this system. The first reason is that the dy- 
namics of spin-1 spinor BECs induced by counterflow 
exhibit characteristic behaviors. Dynamics peculiar to 
counterflow have been observed in two-component BECs 
[6|, l7j. However, there are distinct differences between 
two-component BECs and spinor BECs. A spinor BEC 
has spin degrees of freedom. The number of particles of 
each component is conserved in a two-component BEC, 
whereas it is not conserved in a spinor BEC because of the 
spin-dependent interaction. Therefore, we expect that a 
spin-1 spinor BEC will exhibit dynamics characteristic of 
not only the counterflow but also the spin degrees of free- 
dom. The second reason for studying counterflow insta- 
bility in spin-1 spinor BECs is that counterflow instability 
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can generate turbulence in spin-1 spinor BECs. As stated 
above, one aim of this study is to investigate the behavior 
of turbulence in spin-1 spinor BECs. The third reason is 
that counterflow in two-components BECs has been ex- 
perimentally investigated |l7lfl8j. We expect that it may 
be possible to experimentally study counterflow of spin-1 
spinor BECs. Hence, we study the counterflow instability 
of spin-1 spinor BECs using the Gross-Pitaevskii (GP) 
equation and the Bogoliubov-de Gennes (BdG) equation. 
We consider counterflow between the m = 1 and m = — 1 
components, where m is the magnetic quantum number. 

Our main purposes are to investigate phenomena char- 
acteristic of the spin degrees of freedom induced by the 
counterflow instability and the properties of a turbulent 
state in spin-1 spinor BECs. For the instability, we in- 
vestigate the pattern of the particle number density and 
the time dependence of the magnitude of the spin density 
vector, etc. For the turbulence of a spin-1 spinor BEC, 
we focus on statistical quantities such as the probabil- 
ity density function (PDF) of the magnitude of the spin 
density vector and the spectrum of the spin-dependent 
interaction energy. 

The results obtained reveal that the behaviors of the 
instability and the turbulent state induced by the coun- 
terflow depend greatly on the sign of the spin-dependent 
interaction. We calculate the dispersion relation ob- 
tained from the BdG equation, the PDF of the magni- 
tude of the spin density vector and the spectrum of the 
spin-dependent interaction energy; these quantities ex- 
hibit different behaviors depending on whether the spin- 
dependent interaction is ferromagnetic or antiferromag- 
netic. As for the spectrum of the spin-dependent inter- 
action energy, the —7/3 power law is clearly found in 
the ferromagnetic case, but not in the antiferromagnetic 
case. The — 7/3 power law can be understood by the 
scaling analysis of the time development equation of the 
spin density vector. 

This paper is organized as follows. Section II describes 
the formulation. In Sec. Ill, we analytically calculate 
the BdG equation. Section IV presents the numerical 
results related to the dynamical instability induced by 
counterflow. The turbulent state in a spin-1 spinor BEC 
is treated in Sec. V. In Sec. VI, we discuss some problems 
of our study. Finally, we summarize the findings in Sec. 
VII. 



II. FORMULATION 



A. Gross— Pitaevskii equation 



d h 2 

ih—^ rn = -—V 2 ^rn-^Con^ rn -\-C 1 ^ S-Smn^n, (1) 

n= — l 

where M is the mass of a particle. The total density n 
and the spin density vector s are respectively given by 



m=— 1 



m,n= — l 

where (Si) mn are the spin-1 matrices: 

[HI • 
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The parameters Co and c\ are the coefficients of the 
spin-independent and dependent interactions for two- 
dimensional system. 

The total energy E is given by 




E = 



E [C("^V 2 )^ m ]dr 



- ndr + - 



s dr. 



(7) 



The spin-dependent interaction energy is the last term 
with the coefficient c\ on the right-hand side of Eq. (7). 
The ground state in a homogeneous system without a 
magnetic field is ferromagnetic for c\ < and polar for 
ci > 0. 

The total particle number and the spin in the z direc- 
tion are conserved in the GP model without a magnetic 
field and the dipole-dipole interaction, namely 



1(N 1 +N + N- 1 )=0, 



(8) 



We consider a spin-1 spinor BEC in a homogeneous 
two-dimensional system at zero temperature because this 
system is easy to study theoretically and is well described 
by macroscopic wave functions ^ m (m = 1,0,-1). For 
simplicity, we do not treat a magnetic field or the dipole- 
dipole interaction. The wave functions ip m then obey the 

GP equation pi III : 



^-AU) = 0, 



(9) 



with N m = J |^ m | 2 <ir. Equations (8) and (9) show that 
Ni has the same time evolution as A^_i and that the 
change of Nq is related to that of N± and iV_i; this is 
important for understanding the instability induced by 
the counterflow (see Sec. IV). 



B. Initial state 

We consider the counterflow between the m = 1 and 
m = — 1 components with a relative velocity Vr = VrC x 
in a homogeneous two-dimensional system, where e x is 
a unit vector along the x direction. The initial state is 
expressed by 



M 0)N 




eX p[i(K VR . r -ft)} 



expH(§V R .r+^)] / 



(10) 



where no is the total density and fi\ and are the 
chemical potentials and are equal to cqUq + MV^/8. We 
use this initial state to investigate the counterflow insta- 
bility and the turbulent state in a spin-1 spinor BEC. 



C. Numerical method 

We use the Crank-Nicholson method to numerically 
calculate the GP equation starting from the initial state 
of Eq. (10). The coordinate is normalized by the co- 
herence length £ = h/ \/2Mcqtiq and the box size is 
128 x 128. Space in the x and y directions is dis- 
cretized into 512 x 512 bins. The time is normalized 
by r = H/coriQ. We add some small white noise to the 
initial state of Eq. (10); without this noise, the instabil- 
ity cannot be generated. In our calculations, the noise 
forms the particles of the m = component, which ac- 
counts for 0.1 ~ 0.3% of the total particle number. This 
is consistent with experimental results (2lj . 



III. BDG EQUATION AND DISPERSION 
RELATION 

In this section, we consider a small deviation 5ijj m from 
the initial state of Eq. (10), whose dispersion relation can 
be obtained by linear analysis. 

We can write the wave functions as 



(11) 



Our system is homogeneous so that we can express the 
small deviation 5ip m by plane waves as 



Sifj m (u, 



i(k-r-ujt) _ * -i(k-r-ujt) 

1 ti (7™ C 



)e 



-iA n 




ft-mV R .r 



(12) 



(13) 



where fio is (/ii + /i_i)/2. Substituting Eqs. (11)— (13) 
into Eq. (1) and neglecting the quadratic terms of the 
small deviation, we obtain the following equations: 
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FIG. 1: Imaginary part of the dispersion relation for the an- 
tiferromagnetic case with co/a = 20 and Co > 0. The upper 
and lower graphs are respectively the imaginary parts of cji,-i 
and ujq. The vertical and horizontal axes are respectively the 
relative velocity and wave number in the x direction, which 
is normalized by the sound velocity c s = \J cono/2M and the 
coherence length £ = h/ \/2Mcono. The imaginary part of cjo 
has finite values for any finite relative velocity Vr, while that 
of cji,_i has finite values only for a relative velocity larger 
than the critical value V c . 
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FIG. 2: Imaginary part of the dispersion relation for the ferro- 
magnetic case with co/a = —20 and Co > 0. The upper and 
lower graphs are respectively the imaginary parts of cji,-i 
and ujq. The vertical and horizontal axes are respectively the 
relative velocity and wave number of the x direction, which 
is normalized by the sound velocity c s = \J cqtiq/2M and the 
healing length £ = h/ \/2Mcono. The imaginary parts of uji,-i 
and ujq have finite values for any relative velocity Vr. There- 
fore, the initial state is unstable even for Vr — 0. This is in 
contrast to Fig. 1, which shows the antiferromagnetic case. 



with €k 
fe/2. 



h 2 k 2 /2M and h± = e k +no(co+c 1 )/2±hV R • 



It follows that the small deviations of the m = ±1 com- 
ponents couple with each other, but that the small devi- 
ation of the m = component develops independently of 
the small deviation of the m = ±1 components. This is 
because the initial state of Eq. (10) does not have any 
m = component. Since Eqs. (14) and (15) are eigen- 
value problems, we can obtain the dispersion relations: 



(ficji _i) 2 = e 2 + (c + ci)n e/e + -(Vr • hk) 2 



2\2 



2^2 



1^0 



(18) 



±y (Vr • Hk) 2 e k (e k + c n + cin ) + ra§(c - Ci) 2 e|{19) 

where cjo and o;i 5 _i are the eigenfrequencies of Eq. (14) 
and (15), respectively. 

The dispersion relations of Eqs. (18) and (19) have 
dynamically unstable regions where the imaginary parts 
of cji 5 _i and ujq become finite. Figures 1 and 2 show 
the imaginary parts of the dispersion relations for the 
antiferromagnetic (co/ci = 20, Co > 0) and ferromag- 
netic cases (cq/ci = — 20, Co > 0), respectively. Here, 
the wave number and the velocity are normalized by the 
coherence length £ = h/ \/2Mcqtiq and the sound veloc- 
ity c s = ^Jc^n^j2M . Experiments typically use 23 Na 
and 87 Rb atoms. The interaction parameters for 23 Na 
atoms satisfy cq/ci ~ 20 and cq > 0, and those for 87 Rb 
atoms satisfy cq/ci ~ —200 and cq > 0. If we had used 
~ —200 and Co > 0, it would take much longer 
for the instability to occur. In this study, to extract the 
dynamics characteristic of the ferromagnetic interaction, 
we use cq/ci = —20 and Co > 0. 

In the following, we explain the character of the dis- 
persion relations and show the kinds of dynamics that 
they are expected to give, which can be confirmed by nu- 
merical calculations based on the GP equation (see Sec. 
IV). 

In the antiferromagnetic case, the imaginary part of 
ujo has finite values for all relative velocities Vr except 
Vr = 0, whereas the imaginary part of uoi^-i has finite 
values only for relative velocities larger than some critical 
value V c . The critical velocity V c /c s of in Fig. 1 

(i.e., the lowest relative velocity for which the imaginary 
part of uo\ : -i is finite) is 2^2c\j c§ ~ 0.63. For < Vr < 
V c , only ujo has an imaginary part. This means that the 
density is modulated only for the m = component. 
However, we expect that the density will be modulated 
for the m = ±1 components too. The instability occurs 
for the m = ±1 components even though the imaginary 
part of cji } _i in Fig. 1 vanishes because Eqs. (8) and (9) 
show that increasing the particle number of the m — 
component reduces the particle number of the m = ±1 
components. Consequently, the instability of the m = ±1 
components can occur. When the relative velocity Vr 
exceeds V c , the amplitude of the imaginary part of uji^-i 
is larger than that of cjo- Thus, the instability of the 
m = ±1 components occurs faster than that of the m = 
component. 

In the ferromagnetic case, the dispersion relations for 
oji^-i and ujq have finite imaginary parts at any arbitrary 
relative velocity Vr. This is in contrast to the antiferro- 
magnetic case and it implies that the initial state of Eq. 
(10) is unstable even without a counterflow. In the case 
Vr ~ 0, the instability of all components occurs nearly 
at the same time because the amplitude of the imaginary 
parts of cji ? _i is almost same as that of ujq. On the other 
hand, in the case Vr > 0, the instability of the m = ±1 
components grows faster than that of the m = compo- 
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FIG. 4: Time dependence of the quantities < Sf > (i = 
t,x,y,z) and the particle number Ni (i = t,x,y, z) for the an- 
tiferromagnetic case. Sf (i — t,x,y,z) is defined by Eqs.(20) 
and (21). These results are for numerical calculations with 
co/ci = 20, c > 0, and V R /c s = 0.39. 



IV. COUNTERFLOW INSTABILITY 

We investigate the dynamics of a spin-1 spinor BEC 
induced by counterflow by performing numerical calcula- 
tions based on the GP equation. The dynamics is mainly 
classified according to whether the spin-dependent inter- 
action is antiferromagnetic or ferromagnetic. In this sec- 
tion, we present the detailed dynamics for both cases. 



A. Antiferromagnetic interaction case 



FIG. 3: Density profiles of the m — 1,0,-1 components for 
the antiferromagnetic case at t/r = (a) 0, (b) 315, and (c) 
1500. (d) and (e) profiles of the spin density vector corre- 
sponding to (b) and (c), respectively. The field of view of 
each image is 128£ x 128£. The shading of the arrows in 
(d) and (e) denote the magnitude of the spin density vector. 
These results are for numerical calculations with co/a = 20, 
c > 0, and Vr/c s = 0.39. 



nent, which reflects the amplitude of the imaginary part 
of co>i,_i and ujq. 

Finally, we discuss the isotropy and anisotropy of the 
dispersion relation about the direction of the relative ve- 
locity. Equation (18) shows that uo depends not on the 
direction of the relative velocity but on its amplitude. 
On the other hand, oji^-i depends on the direction of 
the relative velocity as well as its amplitude. Thus, the 
dynamics obtained by the GP equation exhibits differ- 
ent behaviors depending on which eigenfrequencies give 
larger imaginary parts. If ujq (u^-i) is dominant, the 
early dynamics will be isotropic (anisotropic). 

These results obtained from the dispersion relation are 
consistent with numerical calculations based on the GP 
equation (see Sec. IV). 



The dynamics for the antiferromagnetic interaction in- 
duced by counterflow instability are shown. The dynam- 
ics strongly depends on whether < Vr < V c or V c < Vr. 
As shown in Sec. Ill, in the former case, the instabilities 
for all three components are expected to grow simulta- 
neously, whereas the instability for the m = ±1 compo- 
nents is expected to grow faster than that for the m = 
component in the latter case. 

First, we present the dynamics for the case < Vr < 
V c - Figure 3 shows the density profiles of all components 
obtained by numerical calculations with Vr/c s = 0.39. In 
this case, the critical velocity V c /c s is about 0.63. Fig- 
ure 3(a) shows the initial state in which the m = ±1 
components have a velocity Vr relative to the counter- 
flow. Figure 3(b) shows the density profiles at t/r = 315, 
where instabilities occur in all components. The instabil- 
ity in Fig. 3(b) appears as an isotropic density modula- 
tion independent of the direction of the relative velocity 
Vr: the observed density modulation is circular. These 
results can be understood by considering the imaginary 
part of the dispersion relation of Eqs. (18) and (19). For 
< Vr < V c , the instability is induced by the imaginary 
part of the dispersion relation of the m = component, 
which is independent of the direction of the relative ve- 
locity Vr. In addition, as pointed out in Sec. Ill, the in- 
stability of the m = component causes the instability of 
the m = ±1 components through Eq. (8). Thus, the in- 
stability of the m = ±1 components is caused by isotropic 
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FIG. 5: Density profiles of the m — 1,0,-1 components for 
the antiferromagnetic case at (a) t/r = (a) 70 and (b) 1000. 
(c) and (d) Profiles of the spin density vector corresponding 
to (a) and (b), respectively. The field of view of each image is 
128£ x 128£. The shading of the arrows in (c) and (d) denotes 
the magnitude of the spin density vector. These results for 
numerical calculations with co/a = 20, Co > 0, and Vr/c s = 
1.57. 



density modulation of the m = component and thus it 
also exhibits isotropic modulation, as shown in Fig. 3(b). 
Hence, all components exhibit circular density modula- 
tion. After a certain time, the density modulation of all 
components becomes very complicated, as shown in Fig. 
3(c). After Fig. 3(b), the circular modulation expands 
in every component. This results in a narrow path in 
the low density region. Thus, in the transition from Fig. 
3(b) to (c), density modulation with various wave num- 
bers grows with increasing time. As a result, the circular 
density modulation in Fig. 3(b) disappears. 

The behavior of the spin density vector is shown in 
Figs. 3(d) and (e), which correspond to the density pro- 
file in Figs. 3(b) and (c), respectively. The spin density 
vector in Fig. 3(d) almost lies in the x — y plane. In the 
early stages of the instability at t/r = 315, the density 
modulation of the m = 1 component almost overlaps that 
of the m = — 1 component. Thus, the z component of 
the spin density vector is very small. As time progresses, 
the modulation of high wave numbers increases, and the 
m = ±1 components become less overlapped. Therefore, 
the z component grows as shown in Fig. 3(e). 

To understand the time dependence of the magnitude 



of the spin density vector, we numerically calculate the 
following quantities: 



< sf >-- 



dr 



(i = x,y,z), (20) 



<Sf 



>-- 



< sf >, 



(21) 



where A is the area of the system. Their time depen- 
dence is shown in Fig. 4(a). The result of < Sf > means 
that the instability of all components occurs at t/r ~ 300 
because growth of < Sf > and < Sy > requires the insta- 
bility of m = component and that of < Sf > requires 
the instability of the m = ±1 components. Figure 4(b) 
shows that the particle number of the m = component 
increases rapidly at t/r ~ 300, which corresponds to the 
occurrence of density modulation in Fig. 3(b). Hence 
the instability induced by the counterflow starts to ex- 
change the particle number among the three components 
and causes the spin density vector to increase. 

Next, we present the dynamics for V c < Vr. Figure 
5 shows the density profile for Vr/c s = 1.57. Figures 
5(a) shows that the density modulation of the m = ±1 
components is much greater than that of the m = com- 
ponent because Im[u;i } _i] is larger than Im[cJo] m Fig. 1. 
This density profile differs from that for < Vr < V c be- 
cause the density modulation in Fig. 5(a) is anisotropic. 
This is understood by the dispersion relation for uji^-i in 
Eq. (19) which depends on the direction of the relative 
velocity Vr. The low density region of the m = ±1 com- 
ponents in Fig. 5(a) is nucleated due to the growth of 
the density stripe perpendicular to the relative velocity 
Vr. The interval of the stripe is consistent with the most 
unstable wave number obtained for Im[o;i ? _i] in Fig. 1. 
Through this stripe in the low density region in Fig. 5(a), 
the phase of each wave function rapidly changes by about 
7r, whose structure is similar to solitons in one-component 
BECs. This soliton-like structure soon collapses and the 
density modulation of the m = ±1 components become 
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FIG. 6: Time dependence of quantities < Sf > (i = £, x, y, z) 
and particle number N% (i = t,x,y,z) for the antiferromag- 
netic case. Sf (i — t,x,y,z) is defined by Eqs. (20) and 
(21). These results for numerical calculations with co /ci =20, 
Co > 0, and Vr/c s — 1.57. 
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complicated, as shown in Fig. 5(b). Through this col- 
lapse, the density modulation of the m = ±1 components 
becomes isotropic. Similar dynamics to that shown in 
Figs. 5(a) and (b) has been reported in two-component 
BECs p, 0, [13, [HI • Even after a long time, the density 
modulation of the m = component does not increase. 
However, this result depends on the initial noise (see Sec. 
IV C). 

The spin density vector behaves very differently from 
the case when < Vr < V c , as shown in Figs. 5(c) and 
(d). In Fig. 5(c), density modulation occurs only in the 
m = ±1 components (Fig. 5 (a)), so that the spin density 
vector is oriented in only ±z directions. In addition, 
the vector in Fig. 5(c) exhibits the stripe structure as 
that in the density profiles in Fig. 5(a). As time passes, 
the structure in Fig. 5(c) collapses because the density 
profiles lose the stripe structure as shown in Fig. 5(b). 
Then, as the density modulation of the m = component 
does not increase, the spin density vector cannot lean in 
the x — y plane (Fig. 5(d)). Therefore, domains in which 
the spin density vector points in the +z or — z direction 
are formed through the instability and they move around. 

The behaviors of the magnitude of the spin density 
vector and the particle number of each component for 
V c < Vr differ greatly from those for < Vr < V c too. 
Figure 6 shows the time dependence of < Sf > (i = 
t,x,y,z) and the particle number of each component. < 
Sf > grows rapidly at t/r = 70, but the particle number 
of each component remains almost the same as that of 
the initial state. This means that the instability of the 
m = ±1 components occurs, but that of the m = does 
not, which is consistent with the density profiles in Figs. 
5(a) and (b). This dynamics is almost the same as the 
behavior of the two-component BEC [|, 0, E3, EU - 
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FIG. 7: Density profiles of the m — 1,0,-1 components in 
the ferromagnetic case at t/r = (a) 175 and (b) 1000. (c) 
and (d) profiles of the spin density vector corresponding to 
(a) and (b), respectively. The field of view of each image is 
128£ x 128£. The shading of the arrows in (c) and (d) indicates 
the magnitude of the spin density vector. These results are 
for numerical calculation with co/a = —200, Co > 0, and 
Vr/cs = 



B. Ferromagnetic interaction case 

This subsection presents the dynamics with the ferro- 
magnetic interaction induced by the counterflow insta- 
bility. In this case, the imaginary parts of lj^-i and 
coo in Fig. 2 do not exhibit a critical velocity, unlike the 
case for the antiferromagnetic interaction. Hence, the dy- 
namics greatly depends on whether the relative velocity 
is Vr ~ or < Vr. As shown in Sec. Ill, in the former 
case, the instability of all components grows, whereas in 
the latter case the instability of the m = ±1 components 
grows faster than that of the m = component. 

We present the dynamics for the case Vr = 0. In 
this case, there is no counterflow, but instability occurs 
because the initial state is unstable even for Vr = 0. 
Figure 7 shows the density profiles of each component 
obtained by numerical calculations with Vr = 0. As 
expected, the density profiles do not exhibit anisotropy 
(Fig. 7(a)) since the initial state with Vr = is isotropic. 
Isotropic density modulation occurs; it corresponds to 
the most unstable wave number in Fig. 2. The instability 
grows considerably, as shown in Fig. 7(b). 



The spin density vector corresponding to the density 
profiles in Figs. 7(a) and (b) are shown in Figs. 7(c) and 
(d), respectively. When the instability occurs at t/r ~ 
175, the density modulation of the m = 1 component 
overlaps with that of the m = — 1 component. Thus, 
the spin density vector almost lies in the x — y plane, as 
shown in Fig. 7(c). With increasing time, the overlap 
decreases, so that the vector points in various directions, 
as shown in Fig. 7(d). This behavior of the spin density 
vector is similar to that for the antiferromagnetic case 
with < Vr < V c (see Figs. 3(d) and (e)). Note that 
the magnitude of the vector in the ferromagnetic case is 
larger than that in the antiferromagnetic case (indicated 
by the shading of the arrows in the vector plots); this is 
discussed later. 

Figure 8 shows the time dependences of < Sf > 
(i = t, x, z) and the particle number. These results 
show that the instabilities for all components occur al- 
most simultaneously. The density profile for each com- 
ponent in Fig. 7 is consistent with these results. 

We present the dynamics for the case < Vr. Fig- 
ures 9 and 10 show numerical results for Vr/c s = 1.96. 
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FIG. 8: Time dependence of quantities Sf (i = t,x,y,z) and 
the particle number Ni (i = t,x,y,z) for the ferromagnetic 
case. < Sf > (i = t,x,y,z) is defined by Eqs. (20) and (21). 
These results are for numerical calculation with co/a = —200, 
co > 0, and Vr/c s = 0. 



Figure 9(a) shows the instability of the m = ±1 com- 
ponents; density modulation of the m = component 
does not occur because Im[co>i 5 _i] > Im[cJo]- The low 
density regions of the m = ±1 components in Fig. 9(a) 
are the soliton-like structure, which collapse (Fig. 9(b)). 
As time progresses, the instability of the m = compo- 
nent develops, as shown in Fig. 9(c). The spin density 
vector for this case is shown in Figs. 9(d) and (e), which 
correspond to Figs. 9(a) and (c), respectively. In the 
early stages of the instability, density modulation occurs 
only in the m = ±1 components, so that the spin den- 
sity vector points in the ±z directions, which is similar 
to Fig. 5(c). As time increases, the m = component 
grows. Thus, the x and y components of the spin den- 
sity vector become large so that the spin density vector 
points in various directions, as shown in Fig. 9(e). These 
behaviors of the vector are consistent with the time de- 
pendences of the < Sf > and the particle number of each 
component shown in Fig. 10. 

These results reveal that there are obvious differences 
in the behaviors of the magnitude of the spin density 
vector for the antiferromagnetic and ferromagnetic cases. 
Figures 4, 6, 8, and 10 show that in the antiferromagnetic 
case < Sf > tends to decrease, whereas it tends to in- 
crease in the ferromagnetic case. This gives rise to the 
different behaviors of the PDF of the magnitude of the 
spin density vector, which is very important for the tur- 
bulence of a spin-1 spinor BEC. The details are described 
in Sec. V. 



C. Dependence of the dynamics on the initial noise 

We discuss the dependence of the dynamics on the 
initial noise. In our numerical calculations, small white 
noise is added to the initial state. Different samples of 
the white noise are available, although the magnitude is 
fixed. Our numerical results can be classified into four 
categories: (I) antiferromagnetic case with a small rela- 
tive velocity (Fig. 3); (II) antiferromagnetic case with a 
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FIG. 9: Density profiles of the m — 1,0,-1 components for 
the ferromagnetic case at t/r = (a) 40, (b) 125, and (c) 500. 
(d) and (e) profiles of the spin density vector corresponding 
to (a) and (c), respectively. The field of view of each image is 
128£ x 128£. The shaping of the arrows in (d) and (e) indicates 
the magnitude of the spin density vector. These results are 
for numerical calculations with co/a — —200, Co > 0, and 
Vr/cs = 1.96. 



large relative velocity (Fig. 5); (III) ferromagnetic case 
with a small relative velocity (Fig. 7); (IV) ferromagnetic 
case with a large relative velocity (Fig. 9). We performed 
some numerical calculations and we empirically observed 
a strong dependence on the noise sample only for (II). 
This strong dependence is related to the growth of the 
m = component. In this case, the instability of the 
m = ±1 components occurs first (Figs. 5 and 6), which 
is independent of the noise. However, numerical calcula- 
tions reveal that as time increases, the m = component 
may grow depending on the noise sample used. We can- 
not control whether growth occurs or not. In this paper, 
we consider the results that do not depend on the initial 
noise; thus, Figs. 5 and 6 show only the results obtained 
before the m = component grows. 

For cases (I), (III), and (IV), the dynamics is quali- 
tatively unchanged even when the initial noise sample is 
varied. 
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FIG. 10: Time dependences of the quantities < Sf > (i — 
t,x,y,z) and the particle number Ni (i — t,x,y,z) for the 
ferromagnetic case. Sf (i — t,x,y,z) is defined by Eqs. (20) 
and (21). These results are for numerical calculations with 
c / Cl = -200, c > 0, and Vr/c s = 1.96. 



V. TURBULENCE IN A SPIN-1 SPINOR BEC 

We find that, in the ferromagnetic case, the spectrum 
of the spin-dependent interaction energy in the turbulent 
state obeys the —7/3 power law, whereas that in the an- 
tiferromagnetic case does not. This result was obtained 
by numerical calculations and it can be understood in 
terms of scaling analysis. This section mainly considers 
the spectrum of the spin-dependent interaction energy 
and its time dependence. 

The methodology usually used for analyzing classical 
turbulence is applied to a spin-1 spinor BECs in this 
study. Here, we briefly review turbulence in classical flu- 
ids p], Many studies focus on statistical quantities 
and laws because they reflect properties characteristic 
of the complicated motion in turbulence. One quantity 
often investigated is the kinetic energy spectrum. Kol- 
mogorov proposed that the kinetic energy spectrum of 
incompressible fluid obeys a —5/3 power law in fully de- 
veloped homogeneous isotropic turbulence; this can be 
demonstrated by making several assumptions. One as- 
sumption is that the kinetic energy flux is independent 
of the wave number. This means that the kinetic energy 
is transported to a high wave number with a constant en- 
ergy flux in the wave number region that obeys the —5/3 
power law. This result has been confirmed by many nu- 
merical calculations and experiments. The inert ial term 
in the Navier-Stokes equation was found to play a domi- 
nant role in energy transfer. These results show the self- 
similarity of the velocity field in wave number space. On 
the other hand, in real space, a Richardson cascade is 
believed to occur in which large vortices become smaller 
through reconnect ions of vortices. However, this has not 
been confirmed yet. 

We focus on the spectrum of the spin-dependent inter- 
action energy in the turbulence of a spin-1 spinor BEC. 
In analogy with classical turbulence, the flux of the spin- 
dependent interaction energy is expected to be indepen- 
dent of the wave number, which gives rise to the spectrum 
characteristics of this system. 



We derive an expression for the spectrum of the spin- 
dependent interaction energy. The spin-dependent inter- 
action energy E s per unit area is given by 

Es = IaJ s{r)2dr - (22) 

We expand the spin density vector s(r) with plane waves: 
s(r) = J2me ik - r . (23) 



The spin-dependent interaction energy E s is represented 
by s(k) as 

ci 
2 



£ s = f 5>~(fc)| 2 . 



(24) 



Therefore, the energy spectrum of spin-dependent inter- 
action energy is given by 

ci 



E s (k) 



2Ak 



/c<|fci|</c+Afc 

where Ak is 2tt/L for a system size L. 

A. Ferromagnetic interaction case 

Our numerical results reveal that the spectrum of the 
spin-dependent interaction energy in the ferromagnetic 
case obeys the power law shown in Fig. 11 (a). This 
spectrum is numerically calculated at t/r = 5000, which 
is a sufficiently long time after the occurrence of the in- 
stability. This is found by the time dependence of < Sf > 
(i = t, x, y, z) shown in Fig. 11 (b). The spectrum in Fig. 
11 (a) has two regions that are separated by the wave 
number 27r/£ s , which corresponds to the spin coherence 
length £ s = h/ ^/2M|ci|no, (i.e., the characteristic scale 
of spin structures such as domain walls and polar core 
vortex [28j). In the low wave number region (k < 27r/£ s ), 
the spectrum obeys the power law, whereas it does not 
in the high wave number region (2ir/^ s < k). The ori- 
gin of this power law in the low wave number region is 
discussed below. 

We will obtain the power exponent of the spectrum of 
the spin-dependent interaction energy in the low wave 
number region (k < 27r/£ s ) using the GP equation. The 
calculation used to obtain Fig. 11 (a) was for the ferro- 
magnetic case, where the magnitude of the spin density 
vector is expected to have some large value in their PDF. 
To confirm this, we calculated the PDF of the magni- 
tude of the spin density vector. The PDF in Fig. 12 
has a sharp peak, indicating that the magnitude of the 
spin density vector tends to be no in the turbulent state. 
This observation allows us to write the macroscopic wave 
function in the turbulent state with the ferromagnetic in- 
teraction as 




i(0- 7 ) 



cos 2 f N 



^sin/3 
e ia sin 2 § 



(26) 
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and (26): 
d 



&t 



s + (y- V)s 



2M 



s x [V 2 s + (a - V)s], (28) 



2Mn i 



(CVi-iVC) (29) 



with a = (Vno)/no [9l-ll lj . In our case, the total density 
is approximately no in the initial state of Eq. (10), so 
that a vanishes. Therefore, the following equation can 
be used to calculate the power exponent of the spectrum: 



s x V' 2 s, 



(30) 



where space and time are respectively normalized by £ 
and r (t = t/r, V = £V). In the early stages of the 
instability, the wave function has vortices and a soliton- 
like structure. The velocity v can be large in their neigh- 
borhoods. This is similar to the vicinity of vortices in 
one-component BECs, where the velocity is larger than 
the sound velocity in the vortex core region whose size 
is equal to the coherence length. However, the number 
of such defects tends to decrease with time. Therefore, 
the velocity v is lower than the sound velocity c s almost 
everywhere after a sufficiently long time after the insta- 
bility occurs. We confirm that the PDF of the magnitude 
of the velocity has a peak at about l/10th of the sound 
velocity. We thus expect that the term on the right-hand 
side of Eq. (30) is important for transporting the spin- 
dependent interaction energy to a higher wave number. 
This is different from classical turbulence for which the 



FIG. 11: Spectrum of the spin-dependent interaction en- 
ergy E s (k) (upper) and time dependences of the quantities 
< Sf > (i — t,x,y,z) (lower) for the ferromagnetic case at 
t/r = 5000. In the graph of the spectrum, the red squares and 
blue solid and black dashed lines show the numerical results 
for the spectrum, a line proportional to k~ 7 ^ 3 , and the bound- 
ary of the wave number corresponding to the spin coherent 
length, respectively. This spectrum is obtained by numerical 
calculations with cq/ci = —20, cq > 0, Vr/c s = 0.78. 



where a, /?, and 7 are the Euler angles in spin space and 
no is the total density of the initial state. This wave 
function leads to 

s = no (sin j3 cos a, sin j3 sin a, cos j3) . (27) 

In addition, the total density is assumed to be time- 
independent because Co ^> c\. In this condition, the spin- 
independent interaction energy is larger than the kinetic 
and spin-dependent interaction enegies, so that the total 
density n has a weak time dependence. We can obtain 
the time evolution equation of s = s/no from Eqs. (1) 
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FIG. 12: Probability density function of the magnitude of the 
spin density vector for the ferromagnetic case at t/r = 5000. 
This was obtained by performing numerical calculations with 
the following parameters: co/ci = —20, Co > 0, Vr/c s = 0.78. 



11 



(a) 














(b) 


! 




X i 








«■ : \ 




FIG. 13: Time dependence of spectrum of spin-dependent 
interaction energy E s (k) for the ferromagnetic case at t/r = 
(a) 125, (b) 600, (c) 3000, and (d) 5000. In the spectrum, the 
red squares, blue solids, and black dashed lines indicate the 
numerical results for the spectrum, the line proportional to 
A; -7 / 3 , and the boundary of the wave number corresponding to 
the spin coherent length, respectively. This was obtained by 
performing numerical calculations with eo/ci = —20, Co > 0, 
V R /c s = 0.78. 



inertial term in the Navier-Stokes equation is dominant 
for energy transfer. 

To understand the behavior of the spectrum in Fig. 11, 
we consider spin density vector in wave number space. 
We can express Eq. (30) by neglecting the second term 
on the left-hand side by using s = s/uq: 



- kp(fci) x s(k 2 )5 kM+k2 . (31) 

ki,k 2 



In the following, we apply Kolmogorov-type dimensional 
analysis to Eq. (31) [22|, |23| under two assumptions: 
Equation (31) is invariant under the scale transformation 
and the flux of the spin-dependent interaction energy is 
independent of the wave number. We perform the scale 
transformation k — >> £fc, t — >• rjt in Eq. (31). Then, if s is 
transformed to (~ 2 r]~ 1 s^ Eq. (31) will be invariant. We 
can then write the dependence s on k and t as 



s ~ k~ 2 r 



(32) 



We assume that the flux of the spin-dependent interac- 
tion energy is independent of the wave number, which is 
equivalent to assuming the existence of a region in which 
the energy is constantly transported. This result can be 
expressed by 



where e and t s are respectively the energy flux and the 
characteristic time. Using Eqs. (25), (32), and (33), we 
obtain the —7/3 power exponent: 



E s (k) ~ AT^fc-V) 2 



2/s k~ 7/s . 



(34) 



(33) 



This result agrees with the numerical result in Fig. 11, 
where the blue line denotes k~ 7 ^ 3 . The same assumption 
has been used for classical turbulence; the wave number 
region that obeys the —5/3 power law is known as the 
inertial range since it originates from the inertial term in 
the Navier-Stokes equation [1]. However, in our system, 
the inertial term is not important for energy transfer, so 
that the region obeying the power law in Fig. 11 cannot 
be termed the inertial range. 

The power law in the spectrum of the spin-dependent 
interaction energy originates from the fact that the en- 
ergy flux induced by the first term of the right-hand side 
of Eq. (28) is independent of the wave number. This 
nonlinear term contains the second derivative and it dif- 
fers from the first-derivative inertial term in the Navier- 
Stokes equation. The different nonlinear terms in the two 
equations are responsible for the different exponents of 
the energy spectrum. Therefore, the —7/3 power law of 
the spectrum is peculiar to turbulence of a spin-1 spinor 
BEC. 

Figure 13 shows how the —7/3 power spectrum devel- 
ops. In the early stages of the instability, the spectrum 
has a peak corresponding to the most unstable wave num- 
ber of the dispersion relation given by Eqs. (18) and (19), 
as shown in Fig. 13(a). This is confirmed by Fig. 2, from 
which it follows that the most unstable wave numbers kt; 
is approximately equal to 0.3 ~ 0.4. Actually, the den- 
sity modulation at t/r = 125 has a stripe structure that 
resembles those in Figs. 5(a) and 9(a), and the unsta- 
ble wave number corresponds to the wave number of the 
stripe. After the instability occurs, the spectrum changes 
to that in Fig. 13(b). The density then no longer sustains 
the stripe structure and excites modulation of various 
wave numbers. As time increases, the spectrum starts to 
obey the —7/3 power law, as shown in Fig. 13(c). Af- 
ter that, the spectrum continues to obey this power law. 
This is confirmed by Fig. 13(d) showing the spectrum at 
t/r = 5000. 

We calculate the time dependence of the power expo- 
nent by the least-squares method 0, HH • The deviation 
(7 from the straight line obtained by the method is also 
calculated. The results with Vr/c s = 0.39, 0.78, and 
1.96 are shown in Fig. 14, which exhibits that the power 
exponent is approximately —7/3 over a long time. Note 
that since our system is not stationary, the power expo- 
nent asymptotically approaches —7/3 in Fig. 14. Also, 
Fig. 14(a) shows how the time development of the power 
exponent depends on the relative velocity Vr. When the 
relative velocity is small, it takes long time for the value 
of n to approach —7/3. This is confirmed by comparing 
the cases Vr/c s = 0.39 and 0.78. On the other hand, the 
spectrum behaves differently when Vr/c s = 1.96. In this 
case, the value of n approaches —7/3 from above (Fig. 
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14(a)). This can be understood from the imaginary part 
of the dispersion relation in Fig. 2, which shows that 
the unstable region of cji 5 _i is broad for a large relative 
velocity. This means that, unlike the case for a small 
velocity, the instability contains various wave numbers. 
Thus, in this case, the spectrum in the low wave number 
region < k < k s becomes flat in the early stages of the 
instability, leading to the time dependence of n. 

In summary, we observed a —7/3 power law in spin-1 
spinor BECs with the ferromagnetic interaction. This 
law is independent of the relative velocity Vr in our 
numerical calculations based on the GP equation for 
< Vr/c s < 1.96. However, if the relative velocity is 
much greater than the sound speed c s , this law will not 
hold because the total particle density n may be inhomo- 
geneous. 



B. Antiferromagnetic interaction case 

Unlike the ferromagnetic case shown in Fig. 11, the 
spectrum of the spin-dependent interaction energy in the 
antiferromagnetic case does not obey the —7/3 power 
law. This is confirmed by Fig. 15 which shows the spec- 
trum of the spin-dependent interaction energy for the an- 
tiferromagnetic case at t/r = 6000 when the spin density 
vector is greatly disturbed. However, this spectrum may 
show a power law in the narrow range 0.3 < fc£ < 1.4. 
As shown in the following, we cannot estimate the power 
exponent by the simple scaling analysis applied to the 
ferromagnetic case. 

The —7/3 power law for the ferromagnetic case was 
obtained by applying a scaling argument under certain 
assumptions, which are not applicable for the antiferro- 
magnetic case. This can be confirmed by checking the 
PDF of the magnitude of the spin density vector (Fig. 
16). There are two distinct differences between Figs. 12 
and 16: the peak width and the magnitude of the spin 
density vector corresponding to the peak both differ. The 
peak width in Fig. 16 for the antiferromagnetic case is 
larger than that in Fig. 12 for the ferromagnetic case. 
The magnitude of the spin density vector corresponding 
to the peak for the ferromagnetic case is approximately 
no, whereas that for the antiferromagnetic case is smaller 
than no- These results imply that Eq. (26) is not valid 
for the antiferromagnetic case. Thus, Eq. (34), which 
is based on Eq. (26), cannot be applied to the antiferro- 
magnetic case. Therefore, the spectrum for the antiferro- 
magnetic case cannot be analyzed by the simple scaling 
analysis applied to the ferromagnetic case. 



VI. DISCUSSION 

In this section, we discuss the following three topics: 
(i) comparison between turbulence in spin-1 spinor BECs 
and other kinds of turbulence; (ii) possibility of observing 
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FIG. 14: Time dependence of the power exponent and de- 
viation from the straight line obtained by the least-squares 
method for the ferromagnetic case. The value of n in (a) is 
calculated by the least-squares method. The deviation a in 

(b) is defined as yj^2^ =1 (yi — f(xi)) 2 /N, where /, N, and 
(yi,Xi) are the line obtained by the least-squares method, 
the number, and the data set. In this calculation, the wave 
number region (ko < k < k s ) is limited whit ko — 8tt/L. 
This was obtained by performing numerical calculations with 
co/a = -20, c > 0, Vr/c s = 0.39,0.78, 1.96. 



the —7/3 power law, and (iii) how to generate counterflow 
in spin-1 spinor BECs. 



A. Comparison of turbulence in spin-1 spinor 
BECs and other kinds of turbulence 



We compare the turbulence in spin-1 spinor BECs with 
that in one-component BECs and classical incompressible 
fluids. Through this comparison, we discuss some char- 
acteristic properties of turbulence in spin-1 spinor BECs. 
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FIG. 15: Spectrum of spin dependent interaction energy 
E s (k) for the antiferromagnetic case at t/r = 6000. In the 
spectrum, the red squares and blue solid and black dashed 
lines indicate the numerical results for the spectrum, the line 
proportional to k~ 7 ^ 3 , and the boundary of the wave number 
corresponding to the spin coherent length, respectively. This 
was obtained by performing numerical calculations with the 
following parameters: cq/c\ — 20, Co > 0, Vr/c s = 0.78. 



There are some differences between turbulence in one- 
component BECs and that in spin-1 spinor BECs. In the 
former system, the Kolmogorov —5/3 power law spec- 
trum is obtained only if an external force is applied. This 
force injects energy on large scales and generates quan- 
tized vortices. If the external force is weak and the num- 
ber of vortices is not enough, the incompressible kinetic 
energy decreases and the Kolmogorov spectrum tends to 
disappear [26|, HtJ. Thus, the existence of the quantized 
vortices is considered to be important for sustaining the 
Kolmogorov —5/3 power law in the turbulence of one- 
component BECs. On the other hand, the —7/3 power 
law in spin-1 spinor BECs is sustained without the appli- 
cation of an external force, as confirmed by Fig. 14. In 
the ferromagnetic case (see Figs. 8 (a) and 10 (a)), the 
absolute value of the spin-dependent interaction energy 
of Eq. (22) tends to increase without the application of 
an external force due to the ferromagnetic interaction. 
This is in contrast with the reduction in the incompress- 
ible kinetic energy in the turbulence of one-component 
BECs when no external force is applied. We conjecture 
that this increase allows the system to obey the —7/3 
power law even when no external force is applied. 

The self-similarity in wave number space can be 
strongly related to the structures in real space, which is 
very important for understanding the behavior of tur- 
bulence in real space. In the turbulence of the one- 
component BEC, the Kolmogorov —5/3 power law is 
believed to be related to the Richardson cascade of 
quantized vortices, where larger vortices become smaller 



FIG. 16: Probability density function of the magnitude of 
the spin density vector for the antiferromagnetic case at t/r = 
6000. This was obtained by performing numerical calculations 
with the following parameters: co/ci = 20, Co > 0, Vr/c s = 
0.78. 



through reconnections of vortices. This cascade may cor- 
respond to the energy cascade in wave number space, 
where the incompressible kinetic energy is transported 
to a high wave number. We expect that there are some 
spin structures characteristic of turbulence with the —7/3 
power law in a spin-1 spinor BEC. Like the Richardson 
cascade of quantized vortices, some larger spin structures 
may become smaller to transport the energy to a high 
wave number. However, we currently do not know what 
kinds of spin structures are essential for the —7/3 power 
law. 

We compare turbulence in spin-1 spinor BECs with 
that in classical fluids. In classical turbulence where the 
Reynolds number becomes infinite, the inert ial term of 
the Navier-Stokes equation becomes dominant in the in- 
ertial range, which causes the system to obey the —5/3 
power law. We do not know any quantities that corre- 
spond to the Reynolds number in spin-1 spinor BECs. 
However, as our system becomes turbulent, the term 
on the right-hand side of Eq. (30) becomes dominant. 
Hence, the spin-dependent interaction energy in the tur- 
bulence is transported to a high wave number mainly by 
this term; this differs from the turbulence mechanism in 
incompressible classical fluids. This is very significant for 
the —7/3 power law, as pointed out in Sec. V. 

From the above discussion, it follows that the turbu- 
lence in spin-1 spinor BECs has some properties that 
other kinds of turbulence do not have. There are several 
unknowns associated with this turbulence such as the ki- 
netic energy spectrum, the spin structure in real space, 
the interaction between the velocity and the spin field. 
These will be studied in the future. 
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B. Possibility of observing the —7/3 power law 

The possibility of observing the —7/3 power law is dis- 
cussed. We consider that the —7/3 power law may be 
experimentally observed if the following three conditions 
are satisfied. 

The spectrum of trapped systems with the ferromag- 
netic interaction is expected to manifest the —7/3 power 
law less clearly than that the spectrum of homogeneous 
systems. This is because we neglect the term contain- 
ing Vn in Eq. (28) in our derivation of the —7/3 power 
law. Variation of the total particle density n should af- 
fect the spectrum in trapped systems. However, for large 
trapped systems, this is expected to have a small effect 
on the spectrum. Therefore, large trapped systems are 
preferable for observing the —7/3 power law. This is the 
first condition. 

All components of the spin density vector have been 
experimentally observed by a phase contrast imaging 
method [2l|. The expressions for the spectrum of the 
spin-dependent interaction energy show that it is pos- 
sible to obtain a spectrum if the spin density vector is 
observed everywhere. This is because Eq. (25) contains 
only the Fourier component of the vector s(k). Thus, 
the second condition is the observation of the spin den- 
sity vector. 

The experimental resolution in wave number space is 
important for observing the —7/3 power law. Our numer- 
ical calculations reveals that the spectrum of the spin- 
dependent interaction energy obeys the —7/3 power law 
in the low wave number region (k < k s ), as shown in Fig. 
11. In trapped systems, we expect that the —7/3 power 
law holds in the region kR < k < k Sl where kR is the wave 
number corresponding to the system size. Typically, the 
system size is approximately equal to the Thomas-Fermi 
radius. Thus, experiments for observing the —7/3 power 
law need to have a resolution up to the spin coherence 
length £ s in real space. Since Sadler et al. [21] observed 
a polar core vortex, the resolution of their experiments 
seems to be sufficient for observing the —7/3 power law. 
The third condition is the observation with a resolution 
up to the spin coherence length £ s . Note that the larger 
a system size is, the broader the region that obeys the 
— 7/3 power law will be, because kR is inversely propor- 
tional to the system size. Thus, larger systems are more 
suitable for observing the —7/3 power law (this is also 
consistent with the first condition). 

In summary, we conclude that the —7/3 power law may 
be experimentally observed if the above three conditions 
are satisfied. 



C. How to generate count erflow in spin-1 spinor 
BECs 

We discuss the two methods for generating counterflow 
of spin-1 spinor BECs in trapped systems. 

The first method is to use a double well potential. A 



numerical study has already investigated this [H . In this 
case, the m = ±1 components are separately trapped in 
each well. The central barrier between two wells is then 
removed, which generates counterflow in spin-1 spinor 
BECs. 

The second method is to utilize a magnetic field gra- 
dient. This method was used to generate counterflow 
of two-component BECs [6]. In the initial state, the 
m = ±1 components are trapped in a harmonic trap. 
As a magnetic field gradient is applied, one component 
moves in the the field gradient direction, while the other 
component moves in the opposite direction. The mag- 
netic field gradient is switched off when the two compo- 
nents are sufficiently separated, which causes the m = ±1 
components to move toward the center of the trap, gen- 
erating a counterflow. 

In future, we may use these methods to study counter- 
flow in spin-1 spinor BECs in a trapped system. 



VII. CONCLUSIONS 

This paper addressed two main topics: the dynam- 
ics induced by counterflow of the spin-1 spinor BECs in 
a homogeneous two-dimensional system and the turbu- 
lence generated by the counterflow. These themes are 
investigated using the GP and BdG equations. 

The results reveal that the properties of the dynamics 
and turbulence in this system are strongly dependent on 
whether the spin-dependent interaction is ferromagnetic 
or antiferromagnetic. We summarize the results below. 

The dynamics induced by the counterflow in the spin- 
1 spinor BEC was investigated by performing analyti- 
cal calculations using the BdG equation and numerical 
calculations using the GP equation. We obtain the dis- 
persion relations of Eqs. (18) and (19) from the BdG 
equation; these relations show the dynamical instability. 
The dispersion relations depend on the spin-dependent 
interaction, so that Im[cJi 5 _i] and Im[co>o] for the antifer- 
romagnetic case differ from those for the ferromagnetic 
case (Figs. 1 and 2). The numerical calculations reveal 
that, in the early stages of the instability, the dynamics 
can be understood in terms of the dispersion relations. 
The stripe width of the density modulation in Fig. 5(a) 
and 9(a) is approximately equal to the most unstable 
wavelength obtained by the dispersion relations. In ad- 
dition, the isotropy and anisotropy of the density modu- 
lation in Figs. 3(b) and 5(a) can be explained in terms of 
the dispersion relations. The distinct difference between 
the ferromagnetic and antiferromagnetic cases appears 
in the magnitude of the spin density vector, as shown in 
Figs. 3 ~ 10. This is very important for the spectrum of 
the spin-dependent interaction energy in the turbulence, 
as pointed out in Sec. V. These results reveal dynamics 
peculiar to the spin degrees of freedom. 

We studied the turbulence generated by the counter- 
flow in spin-1 spinor BECs by performing numerical cal- 
culations using the GP equation and scaling analysis. 
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In the ferromagnetic case, the spectrum of the spin- 
dependent interaction energy obeys the —7/3 power law 
(Fig. 11). The power exponent —7/3 is obtained by the 
scaling analysis of Eq. (31) when we make the following 
three assumptions: the wave functions are approx- 
imately expressed by Eq. (26), Eq. (31) are invariant 
under the scale transformation, and the flux of the spin- 
dependent interaction energy is independent of the wave 
number. 

On the other hand, for the antiferromagnetic case, the 
spectrum does not exhibit the —7/3 power law. This is 
probably attributable to the breakdown of the validity 
of Eq. (26). In this case, the PDF of the magnitude 
of the spin density vector (Fig. 16) differs from that 
for the ferromagnetic case (Fig. 12), so that such scaling 
analysis used for the ferromagnetic case is not applicable. 



However, in the narrow range of the wave number, the 
spectrum in Fig. 15 may show a power law, which cannot 
be estimated by the scaling analysis. 

These results are characteristic of the spin degrees of 
freedom and are thus not exhibited by turbulence in one- 
component BECs and classical fluids. There are some 
unresolved problems associated with turbulence in spin- 
1 spinor BECs has (see Sec. VI); we intend to study these 
in the near future. 
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